### Number theoretic

import random

def _gcd2(a,b):
	a=abs(a)
	b=abs(b)
	if a < b:
		(a,b)=(b,a)
	while b!=0:
		(a,b)=(b,a % b)
	return a

def gcd(*nums):
	return reduce(_gcd2,nums)

def _lcm2(a,b):
	return abs(a)*abs(b)/_gcd2(a,b)

def lcm(*nums):
	return reduce(_lcm2,nums)

def factorial(n):
	rv=1
	while n>1:
		rv*=n
		n-=1
	return rv

def binomial(n,k): # n!/k!(n-k)!
	if k==0: return 1
	if k>(n+1)//2: return binomial(n,n-k)
	if n==0: return 0
	return n*binomial(n-1,k-1)/k

def multinomial(*coeffs):
	s=0
	rv=1
	for k in coeffs:
		s+=k
		rv*=binomial(s,k)
	return rv

# n^p mod m, for large integers
def pmod(n,p,m):
	n=n%m
	rv=1
	nn=n
	while p>=1:
		if p%2==1:
			rv=rv*nn%m
		p//=2
		nn=nn*nn%m
	return rv
	
# Miller Rabin primality test
# identifies a composite number as prime with probability 4^-k in the choice of the random numbers
def is_prime(n,k=100,rng=random):
	if n==2 or n==3: return True
	if n%2==0: return False
	if n<3: return False
	d=n-1
	s=0
	while d%2==0:
		d/=2
		s+=1
	for _ in xrange(k):
		a=rng.randint(2,n-2)
		x=pmod(a,d,n)
		if x==1 or x==n-1: continue
		f=False
		for r in xrange(1,s):
			x=x*x%n
			if x==1: return False
			if x==n-1: 
				f=True
				break
		if f==False: return False
	return True